<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dlagtf.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="DLAGTF.1"></a><a href="dlagtf.f.html#DLAGTF.1">DLAGTF</a>( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            INFO, N
      DOUBLE PRECISION   LAMBDA, TOL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IN( * )
      DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DLAGTF.19"></a><a href="dlagtf.f.html#DLAGTF.1">DLAGTF</a> factorizes the matrix (T - lambda*I), where T is an n by n
</span><span class="comment">*</span><span class="comment">  tridiagonal matrix and lambda is a scalar, as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     T - lambda*I = PLU,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where P is a permutation matrix, L is a unit lower tridiagonal matrix
</span><span class="comment">*</span><span class="comment">  with at most one non-zero sub-diagonal elements per column and U is
</span><span class="comment">*</span><span class="comment">  an upper triangular matrix with at most two non-zero super-diagonal
</span><span class="comment">*</span><span class="comment">  elements per column.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The factorization is obtained by Gaussian elimination with partial
</span><span class="comment">*</span><span class="comment">  pivoting and implicit row scaling.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The parameter LAMBDA is included in the routine so that <a name="DLAGTF.32"></a><a href="dlagtf.f.html#DLAGTF.1">DLAGTF</a> may
</span><span class="comment">*</span><span class="comment">  be used, in conjunction with <a name="DLAGTS.33"></a><a href="dlagts.f.html#DLAGTS.1">DLAGTS</a>, to obtain eigenvectors of T by
</span><span class="comment">*</span><span class="comment">  inverse iteration.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, A must contain the diagonal elements of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          On exit, A is overwritten by the n diagonal elements of the
</span><span class="comment">*</span><span class="comment">          upper triangular matrix U of the factorization of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LAMBDA  (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          On entry, the scalar lambda.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input/output) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, B must contain the (n-1) super-diagonal elements of
</span><span class="comment">*</span><span class="comment">          T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          On exit, B is overwritten by the (n-1) super-diagonal
</span><span class="comment">*</span><span class="comment">          elements of the matrix U of the factorization of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C       (input/output) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, C must contain the (n-1) sub-diagonal elements of
</span><span class="comment">*</span><span class="comment">          T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          On exit, C is overwritten by the (n-1) sub-diagonal elements
</span><span class="comment">*</span><span class="comment">          of the matrix L of the factorization of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TOL     (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          On entry, a relative tolerance used to indicate whether or
</span><span class="comment">*</span><span class="comment">          not the matrix (T - lambda*I) is nearly singular. TOL should
</span><span class="comment">*</span><span class="comment">          normally be chose as approximately the largest relative error
</span><span class="comment">*</span><span class="comment">          in the elements of T. For example, if the elements of T are
</span><span class="comment">*</span><span class="comment">          correct to about 4 significant figures, then TOL should be
</span><span class="comment">*</span><span class="comment">          set to about 5*10**(-4). If TOL is supplied as less than eps,
</span><span class="comment">*</span><span class="comment">          where eps is the relative machine precision, then the value
</span><span class="comment">*</span><span class="comment">          eps is used in place of TOL.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (output) DOUBLE PRECISION array, dimension (N-2)
</span><span class="comment">*</span><span class="comment">          On exit, D is overwritten by the (n-2) second super-diagonal
</span><span class="comment">*</span><span class="comment">          elements of the matrix U of the factorization of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IN      (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On exit, IN contains details of the permutation matrix P. If
</span><span class="comment">*</span><span class="comment">          an interchange occurred at the kth step of the elimination,
</span><span class="comment">*</span><span class="comment">          then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
</span><span class="comment">*</span><span class="comment">          returns the smallest positive integer j such that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          where norm( A(j) ) denotes the sum of the absolute values of
</span><span class="comment">*</span><span class="comment">          the jth row of the matrix A. If no such j exists then IN(n)
</span><span class="comment">*</span><span class="comment">          is returned as zero. If IN(n) is returned as positive, then a
</span><span class="comment">*</span><span class="comment">          diagonal element of U is small, indicating that
</span><span class="comment">*</span><span class="comment">          (T - lambda*I) is singular or nearly singular,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0   : successful exit
</span><span class="comment">*</span><span class="comment">          .lt. 0: if INFO = -k, the kth argument had an illegal value
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            K
      DOUBLE PRECISION   EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.111"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.112"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="XERBLA.115"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( N.LT.0 ) THEN
         INFO = -1
         CALL <a name="XERBLA.122"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLAGTF.122"></a><a href="dlagtf.f.html#DLAGTF.1">DLAGTF</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      A( 1 ) = A( 1 ) - LAMBDA
      IN( N ) = 0
      IF( N.EQ.1 ) THEN
         IF( A( 1 ).EQ.ZERO )
     $      IN( 1 ) = 1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.137"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
<span class="comment">*</span><span class="comment">
</span>      TL = MAX( TOL, EPS )
      SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
      DO 10 K = 1, N - 1
         A( K+1 ) = A( K+1 ) - LAMBDA
         SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
         IF( K.LT.( N-1 ) )
     $      SCALE2 = SCALE2 + ABS( B( K+1 ) )
         IF( A( K ).EQ.ZERO ) THEN
            PIV1 = ZERO
         ELSE
            PIV1 = ABS( A( K ) ) / SCALE1
         END IF
         IF( C( K ).EQ.ZERO ) THEN
            IN( K ) = 0
            PIV2 = ZERO
            SCALE1 = SCALE2
            IF( K.LT.( N-1 ) )
     $         D( K ) = ZERO
         ELSE
            PIV2 = ABS( C( K ) ) / SCALE2
            IF( PIV2.LE.PIV1 ) THEN
               IN( K ) = 0
               SCALE1 = SCALE2
               C( K ) = C( K ) / A( K )
               A( K+1 ) = A( K+1 ) - C( K )*B( K )
               IF( K.LT.( N-1 ) )
     $            D( K ) = ZERO
            ELSE
               IN( K ) = 1
               MULT = A( K ) / C( K )
               A( K ) = C( K )
               TEMP = A( K+1 )
               A( K+1 ) = B( K ) - MULT*TEMP
               IF( K.LT.( N-1 ) ) THEN
                  D( K ) = B( K+1 )
                  B( K+1 ) = -MULT*D( K )
               END IF
               B( K ) = TEMP
               C( K ) = MULT
            END IF
         END IF
         IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
     $      IN( N ) = K
   10 CONTINUE
      IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
     $   IN( N ) = N
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLAGTF.188"></a><a href="dlagtf.f.html#DLAGTF.1">DLAGTF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
